﻿#!/usr/bin/perl
#modifyId:200712230053
#modify for the key word "origin" go to next line

#modifyId:200712141619
#œ«UÄÉÈë·ÖÎö¡£

#ModifyId:200712012210
# ÓÃ»§¶šÒåÆðÊŒÎ»µã
#

#
#This script is used to create point simulation data file that is mutted
#in every site of the each site.
#
#A T C G
#read the 3utr.fasta file and read the first file
#clean the sequence
#do point muation
#out put the mutation file to pointmut
#Author zuofeng
#Shanghai Center for Bioinformation Technology
#last modification:2007/11/15



use strict;


my $line;
my $title;
my $seq;  # utr sequence

my $infile = shift || die("1 Corret Parameter: pointMut.pl fastaFile start_p end_p outputPath.\n\n -1 for full length");
my $start_p = shift || die("2 Corret Parameter: pointMut.pl fastaFile start_p end_p outputPath.\n\n");
my $end_p = shift || die("3 Corret Parameter: pointMut.pl fastaFile start_p end_p outputPath.\n\n");
my $outpath = shift || die("3 Corret Parameter: pointMut.pl fastaFile start_p end_p outputPath.\n\n#with slash at the end of path"); 
my $out_start_p = $start_p;
if($end_p!=-1)
{
	if( $end_p<0 || $start_p >$end_p)
	{
		die("The start_p should be less than end_p!\n\n");
	}
}


#my $infile = "3utr.fasta";

#read the utr file file in fasta format
unless(open(UTR, $infile))
{
 	die("Read file erro!!!");
}

my $count = 0;
while($line = <UTR>)
{
 chomp($line);
 if($line=~/^>/)
 {   
   if($count<1)
   {
   	$line =~s/>//;
   	$title = $line;
	next;  	 	   			 #continue after get seq name information
   }
   else
   {
   	last; 					 #skip the circle 
   }
   $count++;
 }
 $seq.= $line;
}
$seq = &SeqClean($seq);
close(UTR);

#####################3
print $title;
#
my $startCodon = 0;
my $Cds_Length =0;
my $isStandard = 1;
if($title=~/miGLE\|3\'UTR/)
{
	
	my @field = split(/\|/, $title);
	if($#field == 9)
	{
		print "standard";
		
		$isStandard =1;
		
		$startCodon = $field[7]-1;
		$Cds_Length = $field[9];
	}
	else
	{
		$startCodon = 0;
		$Cds_Length = 0;
		$isStandard = -1;
	}
}
else
{
	$startCodon = 0;
	$Cds_Length = 0;
	$isStandard = -1;
}
####################


if($startCodon <=0 && $Cds_Length <=0)
{
	#ÎªÁËÓëÊý×éµÄÖµ±£³ÖÒ»ÖÂ
	$start_p -=1;
	$end_p -=1;
	
}
else
{
	$start_p = $start_p -($startCodon+ $Cds_Length-1 );
	$end_p  =  $end_p - ($startCodon+ $Cds_Length-1 );
	
	#ÎªÁËÓëÊý×éµÄÖµ±£³ÖÒ»ÖÂ
	$start_p -=2;
	$end_p -=2;
	print $start_p."\n";
	print $end_p;
}
my $outfile = $outpath."pointmut";
unless(open(OUT, ">$outfile"))
{
 die("Create out file erro!!!");
}

#record original sequence
print OUT ">$title|Origin\n";#modifyId:200712230053
print OUT $seq. "\n\n";

my @seqSplit = split(//, $seq);
my $lastPos = $#seqSplit;#»ñµÃÐòÁÐ³€¶È
if($end_p != -1)
{
	if($end_p > $lastPos)
	{
		print ("the end_point is out of limit!\n\n");
		#reset to lastposition
	}
	else
	{
		$lastPos = $end_p;
	}

        if($start_p >$lastPos)#for both very larg situation#200712291331
	{
            $start_p = $lastPos-($end_p-$start_p);
        }
}

if($start_p < 0)
{
	$start_p =0; 
	#³¬œçÖØÖÃ
   if($end_p <0)
   {
       $end_p = 20;
   }
}


my $window = 30;	 # need to confirm;
my $start = 0;
my $end = 0;
my $mutSpace = "ATCG";

for(my $i = $start_p; $i<=$lastPos; $i++)
{
 my $hp = $seqSplit[$i];
 $start = $i - $window;
 $end = $i + $window;
  
 ################set window
 if($start <= 0)
 {
  $start =0;
 }
 #if( $end >= $lastPos)
 #{
  #$end = $lastPos;
 #}

 ############mutation
 if($hp =~/[ATCG]/i)
 {
  #get the sequence 
  my $first = "";
  my $second ="";
  for(my $k = $start; $k<=$end; $k++)
  {
   if($k < $i)
   {
   	$first.=$seqSplit[$k];
   }
   if($k > $i)
   {
   	$second.=$seqSplit[$k];
   }
  } 
  
  #output the mutation data
  my @field = split(//, $mutSpace);
  for(my $m=0; $m<4; $m++)
  {
   print $i." ".$field[$m]."  ".$hp."\n";
   if(uc($field[$m]) eq uc($hp))
   {
   	next;
   }
   if($isStandard >=1)
   {
	my $index = $i + $Cds_Length +1 ;
	print OUT ">".$index."|".$hp."->".$field[$m]."|pointMut\n";   #OUT
	print OUT $first.$field[$m].$second."\n";
   }
   else
   {
	my $index = $i +1;
	print OUT ">".$index."|".$hp."->".$field[$m]."|pointMut\n";   #OUT
	print OUT $first.$field[$m].$second."\n";
   }
   
   
  }
 }
 else
 {
  #This position do not contain atcg, we will skip.
  print "This position is not handled by this program erro/n";
  print "In our program, we only allow point mutation!";
  next;
 }
 print OUT ("\n");
}


close(OUT);

########################################################
sub SeqClean()
{
 my($rawSeq) = @_;

 #for windows sequence
 $rawSeq =~tr/\r//d;

 #If there is space delete all sapce(/d)
 $rawSeq =~tr/ //d; 
 
 #if there is number, delete numbers.
 $rawSeq =~tr/[0-9]//d; 
 
 #modifyId:200712141619
 if($rawSeq=~/[uU]/)
 {
	$rawSeq=~tr/[uU]/T/;
 }
 
 return(lc($rawSeq));
}
